perm filename SAILIB.SAI[SYS,HE]1 blob
sn#004143 filedate 1972-06-08 generic text, type T, neo UTF8
00100 ENTRY INVERSION,BOUNDED,DECOMPOSE,SOLVE,IMPROVE;
00200
00300 BEGIN "SUPPORT"
00400 SAFE INTEGER ARRAY PS[1:50];
00500 SAFE REAL ARRAY R[1:50],DX[1:50];
00600
00700
00800 SIMPLE PROCEDURE SINGULAR(INTEGER WHY);
00900 COMMENT PRINTS ERROR MESSAGES FOR DECOMPOSE AND IMPROVE;
01000 BEGIN
01100 IF (WHY=0) THEN
01200 USERERR(0,1, "MATRIX WITH ZERO ROW IN DECOMPOSE." );
01300 IF (WHY=1) THEN
01400 USERERR(0,1, "SINGULAR MATRIX IN DECOMPOSE. SOLVE WILL DIVIDE BY ZERO." );
01500 IF (WHY=2) THEN
01600 USERERR(0,1, "NO CONVERGENCE IN IMPROVE. MATRIX IS NEARLY SINGULAR." );
01700 END ;
00100 INTERNAL SIMPLE PROCEDURE DECOMPOSE(INTEGER N;SAFE REAL ARRAY A,LU);
00200 COMMENT A,LU[1:N,1:N];
00300 COMMENT USES GLOBAL INTEGER ARRAY PS;
00400 COMMENT COMPUTES TRIANGULAR MATRICES L AND U AND PERMUTATION
00500 MATRIX P SO THAT LU=PA. STORES L-I AND U IN LU.
00600 ARRAY PS CONTAINS PERMUTED ROW INDICES;
00700 COMMENT DECOMPOSE(N,A,A) OVERWRITES A WITH LU;
00800 BEGIN
00900 LABEL ENDKLOOP;
01000 INTEGER I,J,K,PIVOTINDEX;
01100 REAL NORMROW,PIVOT,SIZE,BIGGEST,MULT;
01200 SIMPLE PROCEDURE ILOOP(INTEGER UL;REFERENCE REAL R1,R2);
01300 START_CODE LABEL LP,EU;
01400 MOVE 1,-1('17);
01500 MOVE 2,-2('17);
01600 MOVE 3,-3('17);
01700 SUB 3,K;
01800 JUMPLE 3,EU;
01900 LP: AOJ 1,;
02000 AOJ 2,;
02100 MOVN 4,MULT;
02200 FMPR 4,(1);
02300 FADRM 4,(2);
02400 SOJG 3,LP;
02500 EU: END;
02600
02700 COMMENT INITIALIZE PS,LU AND R;
02800 FOR I←1 STEP 1 UNTIL N DO
02900 BEGIN
03000 PS[I]←I;
03100 NORMROW←0;
03200 FOR J←1 STEP 1 UNTIL N DO
03300 BEGIN
03400 LU[I,J]←A[I,J];
03500 IF (NORMROW<ABS(LU[I,J])) THEN NORMROW←ABS(LU[I,J]);
03600 END;
03700 IF (NORMROW≠0) THEN R[I]←1/NORMROW
03800 ELSE BEGIN R[I]←0; SINGULAR(0); END;
03900 END;
04000 COMMENT GAUSSIAN ELIMINATION WITH PARTIAL PIVOTING;
04100 FOR K←1 STEP 1 UNTIL N-1 DO
04200 BEGIN
04300 BIGGEST←0;
04400 FOR I←K STEP 1 UNTIL N DO
04500 BEGIN
04600 SIZE←ABS(LU[PS[I],K])*R[PS[I]];
04700 IF (BIGGEST<SIZE) THEN
04800 BEGIN BIGGEST←SIZE; PIVOTINDEX←I; END;
04900 END;
05000 IF BIGGEST=0 THEN
05100 BEGIN SINGULAR(1); GO TO ENDKLOOP; END;
05200 IF PIVOTINDEX≠K THEN
05300 BEGIN
05400 J←PS[K];PS[K]←PS[PIVOTINDEX];PS[PIVOTINDEX]←J;
05500 END;
05600 PIVOT←LU[PS[K],K];
05700 FOR I←K+1 STEP 1 UNTIL N DO
05800 BEGIN
05900 LU[PS[I],K]←MULT←(LU[PS[I],K]/PIVOT);
06000 IF MULT ≠0 THEN
06100 COMMENT FOR J←K+1 STEP 1 UNTIL N DO
06200 LU[PS[I],J]←LU[PS[I],J]-MULT*LU[PS[K],J];
06300 ILOOP(N,LU[PS[I],K],LU[PS[K],K]);
06400 COMMENT INNER LOOP. ONLY COLUMN SUBSCRIPT
06500 VARIES. USE MACHINE CODE IF NECESSARY
06600 FOR EFFICIENCY;
06700 END;
06800 ENDKLOOP:
06900 END;
07000 IF (LU[PS[N],N]=0) THEN SINGULAR(1);
07100 END ;
07200
07300
00100 INTERNAL SIMPLE PROCEDURE SOLVE(INTEGER N;SAFE REAL ARRAY LU,B,X);
00200 COMMENT LU[1:N,1:N],B,X[1:N];
00300 COMMENT USES GLOBAL SAFE INTEGER ARRAY PS;
00400 COMMENT SOLVES AX=B USING LU FROM DECOMPOSE;
00500 BEGIN
00600 INTEGER I,J;
00700 REAL DOT;
00800 SIMPLE PROCEDURE ILOOP(INTEGER LL,UL;REFERENCE REAL R1,R2);
00900 START_CODE LABEL LP,EU;
01000 MOVE 1,-1('17);
01100 MOVE 2,-2('17);
01200 MOVE 3,-3('17);
01300 SUB 3,-4('17);
01400 SETZ 4,;
01500 JUMPL 3,EU;
01600 LP: MOVE 5,(1);
01700 FMPR 5,(2);
01800 FADR 4,5;
01900 AOJ 1,;
02000 AOJ 2,;
02100 SOJGE 3,LP;
02200 EU: MOVEM 4,DOT;
02300 END;
02400
02500 FOR I←1 STEP 1 UNTIL N DO
02600 BEGIN
02700 COMMENT DOT←0;
02800 COMMENT FOR J←1 STEP 1 UNTIL I-1 DO
02900 DOT←DOT+LU[PS[I],J]*X[J];
03000 ILOOP(1,I-1,LU[PS[I],1],X[1]);
03100 X[I]←B[PS[I]]-DOT;
03200 END;
03300 FOR I←N STEP -1 UNTIL 1 DO
03400 BEGIN
03500 COMMENT DOT←0;
03600 COMMENT FOR J←I+1 STEP 1 UNTIL N DO
03700 DOT←DOT+LU[PS[I],J]*X[J];
03800 ILOOP(I+1,N,LU[PS[I],I+1],X[I+1]);
03900 X[I]←(X[I]-DOT)/LU[PS[I],I];
04000 END;
04100 COMMENT AS IN DECOMPOSE, THE INNER LOOPS INVOLVE ONLY THE COLUMN
04200 SUBSCRIPT OF LU AND MAY BE MACHINE CODED FOR EFFICIENCY;
04300 END;
04400
04500
04600
00100 INTERNAL SIMPLE PROCEDURE IMPROVE(INTEGER N;SAFE REAL ARRAY A,LU,B,X;REFERENCE REAL DIGITS);
00200 COMMENT A,LU[1:N,1:N],B,X[1:N];
00300 COMMENT A IS THE ORIGINAL MATRIX, LU IS FROM DECOMPOSE,B IS THE
00400 RIGHT HAND SIDE,X IS THE SOLUTION FROM SOLVE. IMPROVES
00500 X TO MACHINE ACCURACY AND SETS DIGITS TO THE NUMBER
00600 OF DIGITS OF X WHICH DO NOT CHANGE;
00700 COMMENT MACHINE DEPENDENT QUANTITIES INDICATED BY 0-0;
00800 BEGIN
00900 LABEL CONVERGED;
01000 INTEGER ITER, ITMAX,I;
01100 REAL RT,T,NORMX,NORMDX,EPS;
01200 EXTERNAL REAL SIMPLE PROCEDURE ADPFOR(INTEGER N;REAL ARRAY A;INTEGER I;REAL ARRAY X;REAL EXTRATERM);
01300 EPS←1.0@-8;
01400 ITMAX←16;
01500 NORMX←0;
01600 FOR I←1 STEP 1 UNTIL N DO
01700 IF (NORMX<ABS(X[I])) THEN NORMX←ABS(X[I]);
01800 IF NORMX=0 THEN GO TO CONVERGED;
01900 FOR ITER ←1 STEP 1 UNTIL ITMAX DO
02000 BEGIN
02100 FOR I←1 STEP 1 UNTIL N DO
02200 R[I]←ADPFOR(N,A,I,X,B[I]);
02300 SOLVE(N,LU,R,DX);
02400 NORMDX←0;
02500 FOR I←1 STEP 1 UNTIL N DO
02600 BEGIN
02700 T←X[I];
02800 X[I]←X[I]+DX[I];
02900 IF (NORMDX<ABS(X[I]-T)) THEN NORMDX←ABS(X[I]-T);
03000 END;
03100 IF (NORMDX≤EPS*NORMX) THEN GO TO CONVERGED;
03200 END ;
03300 COMMENT ITERATION DID NOT CONVERGE;
03400 SINGULAR(2);
03500 CONVERGED:
03600 END;
00100 INTERNAL SIMPLE PROCEDURE INVERSION(REAL ARRAY RES,MAT);
00200 BEGIN "INVERT"
00300 SAFE OWN REAL ARRAY LU[1:4,1:4],B,X[1:4];
00400 INTEGER I,J;
00500 REAL DIGITS;
00600 DECOMPOSE(4,MAT,LU);
00700 FOR J←1 STEP 1 UNTIL 4 DO BEGIN
00800 FOR I← 1 STEP 1 UNTIL 4 DO B[I]←IF I=J THEN 1.0 ELSE 0.0;
00900 SOLVE(4,LU,B,X);
01000 IMPROVE(4,MAT,LU,B,X,DIGITS);
01100 FOR I←1 STEP 1 UNTIL 4 DO RES[I,J]←X[I] END;
01200 END "INVERT";
00100 INTERNAL BOOLEAN SIMPLE PROCEDURE BOUNDED(REAL X,Y;SAFE REAL ARRAY P;INTEGER N);
00200 BEGIN INTEGER I,C,N2,II;
00300 REAL M1,M2,RV,RH;
00400 LABEL NEXTL,ADD4;
00500 C←0;
00600 FOR I←N STEP 2 UNTIL 2*(N-1) DO
00700 BEGIN RV←(P[(I+1) MOD N]-Y)*(Y-P[(I+3) MOD N]);
00800 IF RV<0.0 THEN GO TO NEXTL;
00900 IF RV> 0.0 THEN
01000 BEGIN RH←(P[I MOD N]-X)*(X-P[(I+2) MOD N]);
01100 IF RH < 0.0 THEN
01200 BEGIN C←IF X<P[I MOD N] THEN C+1 ELSE C;
01300 GO TO NEXTL
01400 END;
01500 IF RH>0.0 THEN
01600 BEGIN M1←(P[(I+1) MOD N]-P[(I+3) MOD N])/
01700 (P[I MOD N]-P[(I+2) MOD N]);
01800 M2←(IF P[(I+1) MOD N]>P[(I+3) MOD N] THEN
01900 (P[(I+1) MOD N]-Y)/(P[I MOD N]-X) ELSE
02000 (P[(I+3) MOD N]-Y)/(P[(I+2) MOD N]-X)) ;
02100 IF M1=M2 THEN RETURN (TRUE);
02200 C←IF M1>M2 THEN C+1 ELSE C;
02300 GO TO NEXTL
02400 END;
02500 IF P[I MOD N]=P[(I+2) MOD N] THEN RETURN (TRUE);
02600 IF P[I MOD N]=X THEN
02700 BEGIN C←IF X<P[(I+2) MOD N] THEN C+1 ELSE C
02800 END ELSE C←IF X<P[I MOD N] THEN C+1 ELSE C;
02900 GO TO NEXTL
03000 END;
03100 IF P[(I+1) MOD N]=P[(I+3) MOD N] THEN BEGIN IF (P[I MOD N]-X)*(X-P[(I+2) MOD N])≥0.0
03200 THEN RETURN (TRUE) ELSE GO TO NEXTL END;
03300 II← IF Y=P[(I+1) MOD N] THEN I ELSE I+2;
03400 IF X>P[II MOD N] THEN GO TO ADD4;
03500 IF X<P[II MOD N] THEN
03600 BEGIN C←IF (P[(II+3) MOD N]-P[(II+1) MOD N])*
03700 (P[(II+1) MOD N]-P[(II-1) MOD N])>0.0 THEN C+1 ELSE C;
03800 GO TO ADD4
03900 END;
04000 RETURN (TRUE);
04100 ADD4: I←I+2;
04200 NEXTL: END;
04300 RETURN(IF C MOD 2 =1 THEN TRUE ELSE FALSE)
04400 END;
04500 END "SUPPORT"